import numpy as np

LABEL_DATA_PATH = "/disk527/DataDisk/a804_cbf/datasets/lunar_crater_database_robbins_2018_bundle/data/lunar_crater_database_robbins_2018.csv"

dict = {}
with open(LABEL_DATA_PATH, "r") as f:
    f.readline()
    for line in f:
        diameter = line.split(",")[5]
        diameter = int(float(diameter))
        if diameter in dict:
            dict[diameter] += 1
        else:
            dict[diameter] = 1

dia_list = []
with open("diameter_statistic.csv", "w") as f:
    for key in dict:
        f.write(f"{key},{dict[key]}\n")
        dia_list.append((key, dict[key]))
dia_list = np.array(dia_list)
dia_list = dia_list[dia_list[:, 0].argsort()]
dia_res_cumsum = dia_list.sum() - np.cumsum(dia_list[:, 1])
dia_res_cumsum = np.hstack((dia_list[:, 0, None], dia_res_cumsum[:, None]))


def expect_N(r):
    return dia_list[dia_list[:, 0] >= r, 1].sum()


def radius(theta):
    """theta: 弧度"""
    R = 1737.4
    h = 300
    return np.sqrt(
        h * h + 2 * h * R + 2 * R * R * np.sin(theta) - 2 * R * h * np.cos(theta)
    )


# 积分过程：
theta = np.linspace(0, 31, 10000) * np.pi / 180
N = 0
d_theta = theta[1] - theta[0]
for t in theta:
    N += np.sin(t) * d_theta * expect_N(radius(t) / 1000 * 10)
N = N / np.sqrt(3)
print(N)
